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ABSTRACT 

The flow over two different multi-element airfoil configurations is computed using linear eddy vis- 
cosity turbulence models and a nonlinear explicit algebraic stress model. A subset of recently-mea- 
sured transition locations using hot film on a McDonnell Douglas configuration is presented, and the 
effect of transition location on the computed solutions is explored. Deficiencies in wake profile compu- 
tations are found to be attributable in large part to poor boundary layer prediction on the generating 
element, and not necessarily inadequate turbulence modeling in the wake. Using measured transition 
locations for the main element improves the prediction of its boundary layer thickness, skin friction, 
and wake profile shape. However, using measured transition locations on the slat still yields poor slat 
wake predictions. The computation of the slat flow field represents a key roadblock to successful pre- 
dictions of multi-element flows. In general, the nonlinear explicit algebraic stress turbulence model 
gives very similar results to the linear eddy viscosity models. 


1. INTRODUCTION 

1.1 Overview 

The prediction of high-lift (multi-element airfoil) 
flow fields currently represents a difficult challenge for 
the computational fluid dynamics (CFD) and turbulence 
modeling community. Even in two dimensions, state-of- 
the-art CFD codes fail to predict trends with Reynolds 
number or trends with flap/slat rigging changes suffi- 
ciently accurately. 1 Without the capability to consistent- 
ly predict trends using CFD, aircraft designers must rely 
on heuristic techniques and wind-tunnel experiments, 
which themselves present additional difficulties when 
attempting to scale the results up to flight Reynolds 
numbers. ’ 3 ’ 4 

The flow around a multi-element airfoil or wing is in- 
herently complex. Variations in angle of attack and dif- 
ferent slat / flap settings often present very different and 
distinct challenges. For example, for typical landing 
configurations, viscous effects can dominate compress- 
ibility effects near stall, whereas for typical take-off 
configurations compressibility can dominate the stall 
physics. 5 Also, flap separation is often seen at low to 
moderate angles of attack, whereas stall is often caused 
by an unloading of the aft portion of the main element 
due to rapidly spreading and merging shear layers and 
wakes over the flap. 1 

Although high-lift devices work essentially because 
they manipulate the inviscid flow, 2 viscous effects are 
crucial as well. According to Meredith, 6 some of the vis- 
cous features that can affect 2D multi-element systems 
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include: (1) boundary layer transition, (2) shock/bound- 
ary layer interactions, (3) viscous wake interactions, (4) 
confluent wakes and boundary layers, and (5) separated 
flows. Some important insights into the physical pro- 
cesses that govern high-lift aerodynamics were summa- 
rized by A. M. O. Smith in his landmark 1975 paper. 7 In 
particular, he described several effects that contribute to 
improved high-lift characteristics of multiple elements 
over single elements. 

A great body of work on high-lift flows, both experi- 
mental and computational, has built up since that time. 
(See, for example, discussions in the papers by Lynch et 
al, 1 Haines, 2 Woodward and Lean, Thibert et al, 9 Val- 
arezo, 10 and Ying. 5 ) Some experiments have been per- 
formed expressly for the purpose of CFD code 
validation. For example, an AGARD three-element 
take-off configuration, tested in the early 1970’s and 
subsequently used as a part of a battery of AGARD test 
cases, 1 was recently used as a code validation 
challenge 12 by the CFD Society of Canada. Experiments 
have also been performed in the NASA Langley Low 
Turbulence Pressure Tunnel (LTPT) 13,14 on McDonnell 
Douglas three-element landing configurations, used as 
test cases in a CFD Challenge Workshop held at NASA 
Langley in 1993. Several papers have been written that 
describe various aspects of experimental testing on 
these configurations. ' 19 

Since the CFD Challenge Workshop , there have been 
many CFD results reported for the McDonnell Douglas 
three-element configurations, using a variety of numer- 
ical schemes for the Navier-Stokes equations. The in- 
compressible Navier-Stokes code INS2D, using both 
point-matched and overset grids, has been used exten- 
sively to assess the ability of CFD to model trends with 
configuration changes, as well as to assess the effects of 
tunnel walls, grid density, and turbulence model. 20 " 22 It 
has also been coupled with a transition -prediction meth- 
odology. 23 Compressible Navier-Stokes codes have 
been used as well. Jones et al 24 used the CFL3D code 
with overset grids to compute flow over both the 2D 
McDonnell Douglas three-element configurations as 



well as 3D flow over wings with flaps. The structured 
grid codes OVERFLOW and TLNS3D and the unstruc- 
tured-grid codes NSU2D and FXJN2D have also been 
applied to the 2D McDonnell Douglas configura- 
tions. 5,25,26 

1.2 Code Validation Issues 

Overall, CFD comparisons with experiment for the 
2D McDonnell Douglas three-element configurations 
have been somewhat ambiguous. Considering the inher- 
ent complexity of the flow field, CFD results are surpris- 
ingly good, particularly with regard to surface pressures. 
Unfortunately, the trends with Reynolds number and 
configuration changes are not accurate enough to meet 
the needs of wing designers, and the maximum lift coef- 
ficient and the angle of attack at which it occurs tend to 
be overpredicted. However, it is still not clear at this 
time whether the fault is due to (1) inadequate modeling 
in the CFD (turbulence models, transition models), and/ 
or (2) the fact that CFD is not simulating the same con- 
figuration as experiment. 

The former possibility may play a role because turbu- 
lence models vary in the accuracy to which they predict 
free shear flow behavior (such as spreading rate), and 
most in use today are eddy-viscosity one- and two-equa- 
tion models that do not account for curvature effects or 
anisotropy. It is possible that some missing effects such 
as these may be important in certain regions of multi-el- 
ement flow fields. Additionally, most CFD codes pre- 
scribe rather than predict transition locations; these 
locations can have a significant effect on the solu- 
tion. 27,28 Moreover, the transition process itself is usu- 
ally not modeled. Instead, the eddy viscosity is generally 
“switched on” at transition points. More sophisticated 
transition models, such as that employed by Cebeci 29 (in 
an interactive boundary layer code coupled with an alge- 
braic turbulence model), are currently not in wide use. 
The presence of laminar separations bubbles, which are 
usually associated with high peak suctions and steep ad- 
verse pressure gradients, are also generally not ac- 
counted for. 

The latter possibility - that CFD may not be simulat- 
ing the same configurations as experiment - arises par- 
ticularly because of questions regarding deviation from 
two-dimensionality in experimental tests and the appli- 
cability of 2D CFD simulations to such flows. Addition- 
ally, the effects of surface roughness, tunnel turbulence, 
external tunnel devices, and structural model deforma- 
tions need to be considered. Since multi-element airfoil 
flow appears to be very sensitive to all of these parame- 
ters, it is often difficult to determine precisely why a 
CFD simulation is in error from a given experimental re- 
sult. As an example of the experimental sensitivity, a 
three-element airfoil section representing the section at 
59% span of the A3 10 wing was found to give a negative 
trend of C Lm ^ with Reynolds number between 6 mil- 
lion and 16 million for an unpolished model, but a posi- 
tive trend when the model was polished. 30 

1.3 Turbulence Model Studies 

Rogers et al 31 has shown that several different eddy- 
viscosity models tend to be consistent in their prediction 
of the flow over the McDonnell Douglas three-element 
configuration. The differences between the turbulence 


models are generally smaller than the differences be- 
tween computation and experiment. Some further per- 
spective may be gained on the issue of turbulence 
modeling by considering other experimental and com- 
putational results for different configurations. Squire 32 
and Agoropoulos and Squire 33 compared computed ve- 
locity profiles and turbulence quantities with experi- 
mental results over a two-element model designed to 
explore slat-wake / main-wing-boundary-layer interac- 
tion. They describe the merging process by three re- 
gimes: (1) unmerged, when the wake and boundary 
layer are separated by a potential core, (2) initial merg- 
ing, when the outer part of the wake and inner part of the 
boundary layer are unaffected, and (3) full merging re- 
sulting in a new, thicker boundary layer (although full 
merging usually does not occur prior to the main ele- 
ment trailing edge.) Using the incompressible Navier- 
Stokes equations with a k - e turbulence model and an 
algebraic stress model (ASM), they found good agree- 
ment between both models and experiment in the initial 
stages of the merging, with somewhat better results us- 
ing the ASM. Decelerated flows were predicted with the 
same accuracy as those in constant pressure. Squire also 
discussed the fact that stress transport models often pre- 
dict too rapid a relaxation of the flow following distur- 
bances since model constants are generally weighted 
toward equilibrium flows. This results in a too rapid dis- 
appearance of the velocity defect in the wake, as well as 
of the corresponding shear stress peaks. 

In addition to computing the CFD Challenge Work- 
shop landing configurations, Anderson and Bonhaus 34 
computed the flow over a McDonnell Douglas three-el- 
ement take-off configuration with the unstructured code 
FUN2D. They compared with the experimental results 
of Nakayama et al, 35 including turbulent shear stress da- 
ta. Using the one-equation eddy-viscosity turbulence 
model of Spalart and Allmaras, turbulent shear stress 
profiles were found to be in overall reasonable agree- 
ment with the experiment. Peak levels were underpre- 
dicted on the flap, but results in that region may not have 
been sufficiently grid-converged. 

Lien and Leschziner 3 7 and Lien et al 38 computed the 
flow over both a separated single-element ONERA air- 
foil as well as a two-element (main element and flap) 
NLR-7301 airfoil. 39 They compared results using fc-e 
two-equation eddy-viscosity models with full second- 
moment Reynolds stress closure models on the ONERA 
airfoil, and found that the latter models are generally su- 
perior. However, all of the turbulence models performed 
poorly in the far wake, in regions where the turbulence 
production over dissipation is low. The k - e eddy- vis- 
cosity model performed fairly well on the NLR-7301 
airfoil, although wake mixing was incorrect and there 
was insufficient mixing between the wing wake and flap 
boundary layer. This was attributed to the inability of 
the k - e model to represent curvature-turbulence inter- 
action. 

Godin et al 40 also computed the flow over the NLR- 
7301 airfoil, using the one- and two-equation eddy-vis- 
cosity turbulence models of Spalart and Allmaras 36 (S- 
A) and Menter 41 (SST). They found good agreement in 
general with both models, although the SST model tend- 
ed to be more accurate in separated flow regions while 
the S-A model performed better in attached flows and 
wakes, including the merging wake / boundary layer re- 
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gion. Turbulent shear stress profiles appeared to agree 
better in general with experiment than the results of Lien 
et al 38 using the k - e eddy- viscosity model, but velocity 
profiles were similar. Cao and Kusunose, 27 Jasper et 
al, 42 and Fritz 43 also computed this 2D flow, using var- 
ious zero-, one-, and two-equation eddy-viscosity mod- 
els. Generally positive agreement with experimental 
data was obtained for lift, surface pressures, and veloci- 
ty profiles, although the zero-equation (algebraic) mod- 
el used in Jasper et al resulted in poorer predictions of 
boundary layer parameters and stall onset compared 
with a one^equation model. No comparisons with Rey- 
nolds stress data were made in these three references. 

The results for the NLR-7301 airfoil seem to support 
the finding of Rogers et al 31 (for the McDonnell-Dou- 
glas configuration) that, while there are certainly notice- 
able differences at a detailed level, many eddy-viscosity 
turbulence models overall tend to have only relatively 
minor effects on global properties for multi-element air- 
foil predictions. (This conclusion is of course dependent 
on particulars of the flow field in question. For example, 
for flows with boundary layer separation, models such 
as S-A and SST have been shown to be superior to the 
“standard” Jones-Launder version of k - e 4 ,44 ) The 
question of whether more advanced nonlinear turbu- 
lence models would produce significantly better results 
for multi-element airfoils remains open. 

1.4 Focus of this Paper 

Because of the complexity and interdependence of 
many factors in high-lift computing and experimenta- 
tion, it seems inappropriate to focus solely on the subject 
of turbulence modeling without regard for other factors. 
We therefore take a somewhat broader approach. The 
purpose of this paper is to contribute to a greater under- 
standing of issues surrounding CFD validation against 
high-lift wind tunnel tests. As a consequence, it focuses 
both on experimental as well as numerical issues. Spe- 
cifically, the purpose is threefold: (1) introduce a subset 
of transition locations for the 2D McDonnell Douglas 
three-element configuration recently measured in the 
LTPT with hot film, (2) demonstrate the effect of transi- 
tion location on CFD solutions, and (3) explore the ef- 
fect of a nonlinear explicit algebraic Reynolds stress 
turbulence model 45 on the CFD solution for both the 
AGARD take-off configuration 1 1,12 and the McDonnell 
Douglas landing configuration. 15 ' 19 It is hoped that by 
taking this approach, we can not only determine the im- 
pact of modeling the nonlinear Reynolds stress terms, 
but also demonstrate the importance (and difficulty!) of 
accurately modeling the wind tunnel experiment, partic- 
ularly with regard to the transition process. 

2. DESCRIPTION OF THE AIRFOIL 
CONFIGURATIONS AND CFD GRIDS 

Two multi-element configurations are investigated in 
this paper. Both employ one-to-one point connectivity 
across grid zones. Grids with one-to-one point connec- 
tivity insure conservation across boundaries and provide 
improved continuity of grid spacing at zonal interfaces, 
although this type of grid generation is substantially 
more difficult than for overset grids. 

Only fine grids are employed, since extensive grid 
sensitivity studies have been published previously for 
both the CFL3D code 24 as well as for other 


codes 20,2 1,22,26 over similar configurations. These stud- 
ies indicate that velocity profiles are more sensitive to 
grid density than surface pressures and lift coefficient. 
When between 100,000 - 200,000 grid points are used 
for point-matched grids similar to the grids used in the 
current study, further refinement yields only minor dif- 
ferences in the computed velocity profiles. These ear- 
lier studies also indicate the importance of having 
sufficient resolution in the wakes and boundary layers. 
The grids used in the current study have over 100,000 
grid points with minimum normal wall grid spacings 
that yield an average y+ of approximately 0.3 - 0.4 for 
their respective free stream Reynolds numbers. The 
grids are also clustered in the wake regions of each ele- 
ment. 

The first case is a three-element configuration defined 
as case A2 in a battery of AGARD test cases 11 and used 
as a test case for the code validation challenge 12 by the 
CFD Society of Canada in June 1996. Although tested in 
a wind tunnel with solid walls, the AGARD test case is 
performed on a “free air” grid with far field extent of 
about 10 chords. (The original model has a stowed chord 
of c - 0.7635 m, and the tunnel height to chord ratio is 
H/c = 5.19 .) The slat is positioned at an angle of 25°, 
while the single-slotted flap has a moderate deflection 
angle of 20°, typical of a take-off configuration. A close- 
up of the grid used for this case is shown in figure 1 . This 
is a 4-zone grid with 114,908 grid points. The abscissa 
of this figure is given in terms of x/c , and the stowed 
airfoil would have its leading and trailing edges at x/c 
of 0 and 1, respectively. This configuration is computed 
at M - 0.197 and Re = 3.52 x 10 6 based on stowed-ge- 
ometry chord. Transition is tripped on the main element 
at x/c = 0. 1 25 on both the upper and lower surface, and 
is free on the slat and flap. 
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Fig. 1. AGARD configuration grid, with every other grid point 
removed. 

The second configuration is a three-element McDon- 
nell Douglas configuration, 15 ' 19 tested in the LTPT and 
used as a test case in a CFD Challenge Workshop held 
at NASA Langley in 1993. The model has a stowed 
chord of c = 0.5588 m and the tunnel height is 
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H/c - 4.09 . The particular slat and flap settings em- 
ployed are: slat deflection of 30°, slat gap of 2.95%, slat 
overhang of -2.5%, flap deflection of 30°, flap gap of 
1.27%, and flap overhang of 0.25%. This rigging desig- 
nation is referred to as 30P-30N. It is typical of a landing 
configuration. It is computed using both a “free air” grid 
with far field extent of about 15c , as well as a grid that 
models an angle of incidence of 19° with LTPT walls. 
This latter grid extends 15c upstream and 19c down- 
stream of the model to avoid the possibility of inflow / 
outflow boundary influence on the solution, although 
the actual LTPT test section does not extend this far. A 
close-up of this latter grid is shown in figure 2. This is a 
5-zone grid with 138,389 grid points. It is computed at 
M = 0.2 and Re = 9 x l(r based on sto wed-geometry 
chord. In the wind tunnel test, transition is free on each 
of the elements. 
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Fig. 2. McDonnell Douglas 30P-30N grid, with every other grid 
point removed. 

3. DESCRIPTION OF THE TRANSITION TEST 
AND DATA 

A subset of transition data recently taken in the LTPT 
is presented here for the McDonnell Douglas three-ele- 
ment configuration. The hot-film sensors are similar to 
those used by Nakayama et al 46 for a similar configura- 
tion in the same wind tunnel. They are thin nickel films 
on a 0.05mm polyimide substrate. The sensors are ar- 
ranged in a straight streamwise array at 73.4% span, 
whereas most pressure taps are located at 50% span. 
Side-wall suction is used to maintain approximate two- 
dimensionality for the flow. However, this is optimized 
for 16° angle of attack. Generally, up through this angle 
the discrepancies in pressure distribution between the 
50% span row and an auxiliary pressure tap row near the 
hot films are moderate. Above 16°, three-dimensionality 
becomes more pronounced; hence it may be a contribut- 
ing factor in results shown below. Details can be found 
in Bertelrud. 

Three-dimensionality affects the data in two ways. 
First, span-loading variations, which manifest them- 
selves as variations in effective angle of attack, cause 


the locations of suction peaks and stagnation points to be 
displaced in the chord-wise direction. Second, separated 
or low-shear regions, as well as brackets, etc., may cause 
spanwise flow that alters the transition data. Assuming 
experimental data are available, the upper surface 
chordwise effects may be accounted for in the CFD, at 
least in part, by relating the transition to the distance 
downstream of the suction peak as opposed to using an 
absolute geometric location. However, for all the Mc- 
Donnell Douglas cases in this paper, this is not neces- 
sary because the computed suction peak locations agree 
with those in the experiment to within a distance less 
than the densest pitch of the hot films. Spanwise effects 
are not easily accounted for. 

Transition locations are determined based on informa- 
tion from 359 surface hot films on the three elements. 
The start and end of transition are extracted from an 
analysis of a combination of standard deviation, skew- 
ness, and flatness factors of the hot film signals, plus in 
most cases an additional examination of the signal traces 
for auto- or cross-correlations. Note that even though the 
films are located at 2.54mm pitch (or A s/c = 0.0045 ) in 
the most dense regions, the flows at higher angles of at- 
tack have fairly short transition regions, so that only one 
or two films pick up the feature. 


Table 1 Hot film transition data for 30P-30N configuration, 
M=0.2, Re=9 million 



8 deg 

19 deg 

SLAT 

s/c 

x/c 

s/c 

x/c 

lower suet peak 

n/a 


n/a 


lower trans S 

n/a 


n/a 


lower trans E 

n/a 


n/a 


upper suet peak 

.015 

-.084 

.006 

-.0852 L 

upper trans S 

.030 

-.077 

.010 

-.0853 

upper trans E 

.060 

-.057 

.020 

-.082 

MAIN - ‘ 

s/c 

x/c 

s/c 

x/c 

lower suet peak 

-.510 

.526 

-.620 

.635 

lower trans S 

-.405 

.422 

-.670 

.685 

lower trans E 

-.510 

.526 

n/a 


upper suet peak 

.000 

.0496 

.008 

.055 

upper trans S 

.012 

.058 

.016 

.061 

upper trans E 

.025 

.068 

.025 

.068 

FLAP 

s/c 

x/c 

s/c 

x/c 

lower suet peak 

n/a 


n/a 


lower trans S 

n/a 


n/a 


lower trans E 

n/a 


n/a 


upper suet peak 

.026 

.888 

.024 

.886 

(mid) 





upper suet peak 

.020 

.882 

.016 

.878 

(aux) 





upper trans S 

.030 

.892 

.030 

.892 

upper trans E 

.070 

.931 

.060 

.921 


Table 1 lists the starting and ending transition loca- 
tions for the 30P-30N configuration at two angles of at- 
tack, along with locations of the suction peaks. Results 
are given both in terms of the surface coordinate s/c 
(where s/c = 0 corresponds to the location on the for- 
ward part of each element in the stowed position where 
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y/c = 0 ) as well as the deployed x/c coordinate. De- 
pending on location, the data are accurate to between 
A s/c = ±0.002 - 0.005 . Possible film sheet influence 
on transition has not been taken into account. An “n/a” 
for the transition start (S) or end (E) indicates that a de- 
finitive start or end to transition, respectively, was not 
detected prior to the cusp on the lower surfaces of the 
slat and main elements or prior to the trailing edge on the 
flap. The “L” for the upper suction peak at 19° indicates 
that the location is actually on the lower side of the ele- 
ment when the slat is in its rotated position. 

Note the discrepancy present in the location of the flap 
suction peak between the 50%-span station (mid) and 
the auxiliary station near the hot films at 73.4%-span 
(aux) at both angles of attack. Computed suction peaks 
agree best with the “mid” data. This table represents 
only a subset of the transition location measurements 
taken. A complete set of more than 50 cases, including 
different Reynolds numbers, Mach numbers, and flap 
settings as well as more details of the experimental da- 
tabase, can be found in Bertelrud 47 and Bertelrud et al 48 

4. DESCRIPTION OF THE COMPUTER CODE 

The computer code CFL3D 49,50,51 ’ 52 solves the three- 
dimensional time-dependent thin-layer (in each coordi- 
nate direction) Reynolds-averaged Navier-Stokes equa- 
tions with an upwind finite-volume formulation. It can 
solve flows over multiple-zone grids that are connected 
in a one-to-one, patched, or overset manner, and can em- 
ploy grid sequencing, multigrid, and local time stepping 
when accelerating convergence to steady state. Upwind- 
biased spatial differencing is used for the inviscid terms, 
and flux limiting is used to obtain smooth solutions in 
the vicinity of shock waves, when present. Viscous 
terms are centrally differenced. The flux-difference- 
splitting (FDS) method of Roe 53 is employed to obtain 
fluxes at the cell faces. 

The CFL3D code is advanced in time with an implicit 
three-factor approximate factorization method. The im- 
plicit derivatives are written as spatially first-order ac- 
curate, which results in block tridiagonal inversions for 
each sweep. However, for solutions that utilize FDS the 
block tridiagonal inversions are further simplified with 
a diagonal algorithm (with a spectral radius scaling of 
the viscous terms). Further details of the CFL3D com- 
puter code are not given here; they can be found in the 
cited references. CFL3D has been used previously to 
compute 2D multi-element airfoil flow fields with gen- 
erally good success. 24 ’ 42 

When free-air grids are used, a far-field point vortex 
correction 54 is applied at the outer boundary. This cor- 
rection has a significant impact, particularly on the com- 
puted drag for multi-element airfoil flows. 12 When 
internal flow on a grid with walls is computed, the in- 
flow total pressure and total temperature are specified 
according to isentropic flow relations, and the outflow 
back pressure is set to obtain the desired inflow Mach 
number. 


5. TURBULENCE MODELING 

The S-A 36 one-equation and SST 41 two-equation 
eddy-viscosity turbulence models employed in this 
study are documented in their respective references. 
However, the S-A model has a further modification to 
improve convergence by preventing S from going neg- 


ative, suggested by P. Spalart in a private communica- 
tion (March 1993), that 

5 = / v3 s + -n- 2 / v2 0) 

K d 

where 

(2) 

, (I +X/vl)(l-/v2) ,, 3 


and c v2 is taken as 5. The SST model has also been mod- 
ified slightly from its original form. Details are given in 
Menter and Rumsey. (Note that there is an error in 
equation 17 in that reference. It should read:' 

T = /nax(2r 3 ;r,) (4) 

Also, in equation 1 in that reference, the terms <j*p, and 
on, should be replaced by n + o*|i, and n + op,, re- 
spectively.) 

The Gatski-Speziale explicit algebraic stress model 45 
(EASM) represents an effective compromise between 
the full second-moment closure and a two-equation 
eddy-viscosity model. It extracts an algebraic relation- 
ship between the turbulent Reynolds stress and the mean 
velocity field by assuming equilibrium hypotheses on 
both the convective and diffusive terms of the Reynolds 
stress transport equation. The EASM model was first 
implemented in CFL3D by Abid et al 56 in both k-<o 
and k-e form. Further modifications to the model have 
been made since that time. 57 For all the results in this pa- 
per, the model is implemented in k - co formulation as 
follows. 


Dt * p 3jC y L 

K)|] 

(5) 

De> „ a , 2 1 3 

_ = P ._ P<0 

[(-Si 

(6) 

where the production terms are: 


„ M 

P * = - T ^. 

(7) 

P = Y-f-T..— 1 
• k{ vdxj) 

(8) 

and a k = \A^ o M 

v»_fi r ^ Jt ~ 1 ~ ... 

= 2.2, 

_ A A 1 ^ 

P = 0.83 , 

:i:i : 


eddy viscosity term employed in the diffusion terms of 
equations (5) and (6) is given by 


n,* 


c * ^ 
C » (0 


(9) 
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where c M * = 0.081 . 

The explicit nonlinear constitutive equation that is 
used to close the Reynolds-averaged Navier-Stokes 
equations is given (after regularization) by 

( 10 ) 

^^(s ik w tj + s jk w ki ) + ^^{s ik s kj - \s kl s tl S„) 

The turbulent viscosity |x, is 

n, = sf <“) 

and 


c 




3(l+Tl 2 ) + 0.2(Tl 6 -fC 6 ) 
3+ii 2 + 6CV + 6C 2 + il 6 + ; 6 ' 


( 12 ) 


„ ®2/r* o \ 1/2 

^ 


C = §(«W /2 (14) 

and a, = (4/3 - C 2 )(g/ 2) , a 2 = (2 - C 3 )(g/2) , 
a 3 = (2-C 4 )(g/2), and g = (C,/2 + C 5 - 1) . The 
constants governing the Speziale-Sarkar-Gatski (SSG) 
pressure-strain correlation model 58 are: C x = 6.8, 
C 2 = 0.36 , C 3 = 1.25 , C 4 = 0.4 , and C s = 1.88 . 

The p/ terms in equation (10) are given by 


where 


3 (1 ±20 c 
11 3 + T| 2 + 6£ 2 tj 2 + 6£ 2 + T| 6 + 


6, RESULTS AND DISCUSSION 
Results are computed for both the AGARD and the 
McDonnell Douglas configurations. Three turbulence 
models - S-A, SST, and EASM - are used for the 
former, but only S-A and EASM are used for the latter. 


6.1 AGARD Configuration 

For the AGARD configuration, surface pressure coef- 
ficients and boundary layer total pressure coefficient 
profiles are compared with experiment 1 1 at two correct- 
ed angles of attack of 4.01° and 20.18°. The effect of 
transition is not assessed for this configuration. Compu- 
tations are performed with transition set at x/c - 0.125 
on the main element upper and lower surface, and with 
the slat and flap computed “fully turbulent.” The impli- 
cations of these transition settings will be discussed in 
the context of the McDonnell Douglas configuration in 
section 6.2.1. 

Computed surface pressure coefficients are compared 



0.8 0.9 1.0 1.1 12 1.3 


x/c 

(c) flap 

Fig. 3. AGARD surface pressure coefficient, alpha=4<01°. 

with experimental results at an angle of attack of 4.01° 
in figures 3(a) - (c). On the flap and the main element, 
all three turbulence models predict c p for this case in 
good agreement with experimental levels. On the slat, 
all three models predict c p in good agreement with ex- 
periment over the upper surface, and EASM yields sig- 
nificantly different pressures from the other two models 
in the cove region, closer to experiment. 

A similar comparison of surface pressure coefficients 
at a higher angle of attack of 20. 18° is shown in figures 
4(a) - (c). In this case, all turbulence models yield sim- 
ilar results, in good agreement with the experimental re- 
sults. 

Boundary layer profiles of total pressure coefficient 
are given at x/c = 0.35 on the upper surface of the main 
element, as well as at three stations on the upper surface 
of the flap, as shown in figure 5. Results are shown for 
the 4.01° case in figures 6(a) - (d), and for the 20.18° 
case in figures 7(a) - (d). The parameter d in the figures 
represents the normal distance from the airfoil surface. 
At the lower angle of attack, all three models give simi- 
lar results, in good agreement with experiment, with the 
exception that EASM predicts more mixing between the 
flap boundary layer and main element wake than exper- 
iment or the other models. Also, the computations with 
all three turbulence models show evidence of a slat wake 
whereas the experiment does not. It is interesting to note 
that even though EASM shows different slat cove c p 
predictions in figure 3(a), its slat wake does not exhibit 
much difference from the other models. 

At the higher angle of attack, all of the models miss 
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x/c 

(a) slat 


x/c 

(b) main 



Fig. 4. AGARD surface pressure coefficient, alpha=20.18°. 



p to t Ptot 

(a) x/c=0.35 (main) (b) x/c=0.9 1 (flap) 



Fig. 6. AGARD total pressure coefficient profiles, alpha = 4.01°. 


the locations of the main and slat wake at the two aft- 
most stations on the flap. However, as will be shown in 
section 6.2.2, modeling the wind tunnel walls has a large 
effect on the velocity profiles over the flap when the 
multi-element airfoil is at high angle of attack: compu- 
tations with walls included tend to shift the wake loca- 
tions upward. This tendency would improve the current 
predictions. There are also small differences in total 
pressure coefficient levels for the three turbulence mod- 
els there; in particular, EASM yields a deeper slat wake 
deficit than the other two models, while SST yields a 
deeper main-element wake deficit. 



Fig. 5. AGARD profile locations. 

Results for this configuration indicate that, while 
there are minor variations between the three turbulence 
models, all produce very similar surface pressure and 
boundary layer total pressure predictions in general. 
Lift, drag, and moment coefficient over the range of an- 
gle-of-attack through maximum lift are also predicted 



Plot Ptot 

(a) x/c=0.35 (main) (b) x/c=0.91 (flap) 



Fig. 7. AGARD total pressure coefficient profiles, alpha = 
20.18°. 
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Fig. 8. AGARD integrated force coefficients. 

similarly by the three turbulence models for this case, as 
shown in figures 8(a) - (c). The CFD integrated quanti- 
ties compare well with experiment over a large range of 
angles of attack, but computations overpredict the max- 
imum lift coefficient and the angle at which it occurs. 
The computations also do not exhibit slat separation or 
the rapid drop-off of lift past C Lmax typical of take-off 
configurations with leading-edge stall. The drag and 
moment coefficients also do not agree with experiment 
at the highest angles of attack past stall. 

6.2 McDonnell Douglas Configuration 

Most of the following results use the McDonnell Dou- 
glas 30P-30N landing configuration at an angle of attack 
of 19°, although some results are also given for a lower 
angle of attack of 8°. It should be stressed that the hot 
film data and the pressure, velocity, and skin friction 
data were taken in independent wind tunnel tests. 15 " 17,47 
Velocity profiles are given at the locations shown in fig- 
ure 9. When transition locations are set according to the 
experimental data of Table 1, they are assumed to be at 
the ending locations from the hot film data. One excep- 
tion to this is on the slat upper surface at an angle of at- 
tack of 19°. In this case, placing the CFD transition at 
the ending location results in separation, since the lami- 
nar boundary layer cannot negotiate the severe adverse 
pressure gradient. (Although not shown, this result has 
also been confirmed with computations using a finer 
embedded grid around the slat nose, as well as with 
computations using an independent boundary layer 
code. The result is also the same whether or not tunnel 
walls are modeled in the CFD.) Therefore, transition is 
instead placed as far downstream as possible such that 
the laminar boundary layer on the slat does not separate. 



Fig. 9. 30P-30N profile locations. 

6.2.7 Effect of Transition Locations 

The effect of transition locations is explored using the 
S-A turbulence model on the free-air grid only. It is very 
important to remember that the CFL3D code does not 
employ any transition modeling. The code merely ze- 
roes out the production terms in the turbulence model 
equations in the regions where laminar flow is desired. 
There is no transition state at all; the turbulence produc- 
tion is either “off’ or “on.” (Turbulence transport and 
diffusion are always active, however.) This transition 
procedure is different from the “trip source term” proce- 
dure employed by Spalart and Allmaras. The current 
method is used to maintain consistency between the dif- 
ferent turbulence models in CFL3D. 

It is also important to realize that although “fully tur- 
bulent” solutions are accomplished with the production 
terms of the turbulence models active everywhere, this 
does not necessarily guarantee that a turbulent boundary 
layer will exist as far forward as the stagnation point. In- 
stead, in a “fully turbulent” solution, the models essen- 
tially transition on their own. In other words, the 
solution to the turbulence equations results in eddy vis- 
cosity levels which may or may not be large enough to 
create a turbulent boundary layer profile. Generally, at 
high Reynolds numbers the “fully turbulent” computa- 
tions result in transition locations well forward of the lo- 
cations seen experimentally. For example, at 19° the 
locations where the S-A model transitions on its own are 
listed in Table 2. (Transition is assumed to occur when 
the maximum eddy viscosity level nondimensionalized 
by molecular viscosity in the boundary layer first 
exceeds unity.) Stagnation points on the lower surfaces 
of the slat and main elements are x/c = -0.031 and 
0.154 , respectively. The “fully turbulent” transition lo- 
cations are all well forward of the experimental regions 
from Table 1. 


Table 2 “Fully turbulent” transition x/c locations for S-A model 
on 30P-30N configuration, M=0.2, Re=9 million, alpha=19°. 


transition 

x/c, S-A fully turbulent 

x/c range, 
from Table 1 

SLAT lower 

-.034 

n/a - n/a 

SLAT upper 

-.051 L 

-.0853 --.082 

MAIN lower 

.177 

.685 - n/a 

MAIN upper 

.125 L 

.061 - .068 

FLAP lower 
FLAP upper 

turbulent everywhere 
turbulent everywhere 

n/a - n/a 
.892 - .921 
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(a) slat 



x/c 

(b) main 
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(c) flap 


Fig. 10. Effect of transition location on 30P-30N surface pres- 
sure coefficient, S-A model on free-air grid, alpha=19°. 

Figures 10(a) - (c) show surface pressure coefficients 
on each of the elements for both the fully turbulent solu- 
tion as well as a solution for which the transition on each 
of the elements is specified according to the experimen- 
tal data in Table 1, as discussed above. Both computed 
results are in good agreement with experiment, 15,17 with 
only small variations evident on the upper surfaces be- 
tween the fully-turbulent and transition-specified re- 
sults. As pointed out in section 1.1, the multi-element 
interactions largely manipulate the inviscid flow; there- 
fore it is not surprising that transition location has only 
a small effect on the surface pressures. However, these 
small variations in pressure translate to a difference in 
lift coefficient of about 0.12, of the same order of mag- 
nitude as experimental differences in lift coefficient due 
to flap rigging changes that wing designers would like to 
predict with CFD. 1 Also, as will be shown next, the ve- 
locity field is quite sensitive to transition location. Be- 
cause viscous effects are a contributing factor to skin 
friction forces, possible boundary layer separation, and 
the stall process governed by the unloading of the aft 
end of the main element, it is clearly important to predict 
boundary layer and wake profiles accurately. 

To compute the boundary layer and wake thickness 
correctly over the main element, transition locations 
must be set correctly. This effect is demonstrated for the 
upper surface in figures 1 1(a) and (b). The velocity pro- 
files at the stations x/c = 0.45 and x/c = 0.89817 are 


shown for the main element fully turbulent versus main 
element with upper surface transition set close to the 
transition-ending location from Table 1. The slat and 
flap are both fully turbulent in this case. When the main 
element is fully turbulent, its boundary layer profile (be- 
low “slat wake” in figure 11(a)) is too thick compared 
with experiment. 15 When transition is set, the profile be- 
low the slat wake agrees better. The part of the wake due 
to the main element upper surface boundary layer (upper 
part of “main wake” in figure 11(b)) is also better pre- 
dicted with transition specified. In addition, the comput- 
ed skin friction coefficient on the main element, 
although still low at the x/c = 0.45 station, is in better 
agreement with experiment 16 when the transition is set, 
as seen in figure 12. 



(a) x/c=0.45 (main) 



(b) x/c=0.89817 (flap) 

Fig. 11. Effect of main element transition location on 30P-30N 
velocity profiles, S-A model on free-air grid, alpha=19°. 

In figures 1 1(a) and (b), when the slat is fully turbu- 
lent the slat wake is predicted to be too wide and deep. 
Because it is likely that the slat transition location has a 
large influence on the slat wake structure over the down- 
stream elements, we next turn our attention to an inves- 
tigation of its effect. Both the upper and lower surfaces 
of the slat contribute to the shape of the slat wake. How- 
ever, since the lower surface has the added complication 
of separated flow followed by reattachment in the cove, 
we focus our attention solely on the effect of upper sur- 
face transition location. The lower surface transition is 
set at the cusp for all cases shown here. 

Unfortunately, the upper surface transition location 
on the slat cannot be pushed further downstream than 
approximately x/c = -0.0847, while still maintaining 
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x/c 


Fig. 12. Effect of main element transition location on 30P-30N 
upper surface main element skin friction, S- A model on free-air 
grid, alpha=19°. 

an attached laminar boundary layer (experimental tran- 
sition range from Table 1 is x/c = -0.0853 to -0.082 ). 
To effectively force transition further downstream, we 
employ suction over a small region on the slat in the 
computations to prevent the laminar boundary layer 
from separating. Computed velocity profiles are shown 
in figures 13(a) - (c). In these figures, the locations A, 
B, and C correspond with x/c = -0.0847 , 
x/c = -0.082 , and x/c = -0.073 , respectively. Their 
locations relative to the starting and ending locations 
from Table 1 are shown in figure 14. Wake velocity pro- 
files that result from transition at location C agree best 
with experimental profiles, although there are still dif- 
ferences between computation and experiment, particu- 
larly in the slat wake / main boundary layer merging 
region at x/c = 0.85 . In the notation of Squire, 32 (dis- 
cussed in section 1 .3), this latter station is in the region 
of initial merging, whereas the first two stations are in 
the unmerged regime. Note that transition location C is 
downstream of the experimentally-measured transition 
location region. 

Even at the first station x/c = 0.1075 (figure 13(a)), 
which is on the main element just aft of the slat trailing 
edge, the wakes predicted using fully turbulent. A, or B 
have deficits and thicknesses that are too large. (The off- 
set velocity difference between computation and exper- 
iment in this figure is probably a result of improper 
calibration in the experimental data, as suggested by F. 
Spaid in a private communication (May 1997)). Using 
location C, the computed slat wake width and deficit 
roughly agree with experiment at this first station as well 
as over the entire length of the main element. This indi- 
cates that (ignoring the wake / boundary layer merging 
region at x/c = 0.85 for now) the S-A turbulence model 
can do a fairly good job representing the wake develop- 
ment itself given a good initial profile. Therefore, poor 
agreement of the slat wake with experiment may not be 
due to any particular failure of the turbulence model in 
modeling the wake, but rather the fact that the computed 
boundary layer leaving the slat is wrong to start with. 

At a lower angle of attack of 8°, the experiment 15 in- 
dicates an extremely diffuse slat wake at and beyond 



(a) x/c = 0.1075 (main) 



(b) x/c = 0.45 (main) 



(c) x/c = 0.85 (main) 


Fig. 13. Effect of slat transition location on 30P-30N velocity pro- 
files, S-A model on free-air grid, alpha=19°. 

x/c = 0.45 , whereas computations do not. Velocity 
profiles are shown in figures 15(a) and (b) at the first 
two stations of the main element only (results at 
x/c = 0.85 are similar). When the slat is fully turbulent, 
the computed slat wake is too deep and wide at the first 
station x/c = 0.1075, and remains so further down- 
stream. However, even moving the slat transition loca- 
tion to D (x/c = -0.056, which is the approximate 
ending location of transition from the experiment) or to 
E (x/c = 0.015 , near the slat trailing edge) still results 
in too large a wake deficit. (In this case, boundary layer 
suction is not necessary to move the transition location 
downstream, since the adverse pressure gradient is not 
too severe for the laminar boundary layer to handle.) It 
therefore appears that the computations over the slat 
may again be incorrect by the time they leave the slat 
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x/ c 

Fig. 14. 30P-30N upper surface slat transition locations. 


trailing edge. At this angle of attack, however, no upper 
surface slat transition location provides satisfactory slat 
wake profiles in good agreement with experiment. Al- 
though not shown, results are similar using EASM. 

At this time, it is not known what physical processes 
on the slat are not being modeled correctly in the com- 
putations. Although we strongly suspect that the prima- 
ry cause is related to transition (upper and/or lower 
surface transition locations, transition modeling, lami- 
nar bubble modeling, possible transition in the cove 
free-shear layer, etc.), it is also possible that 3D effects 
or possibly unsteady effects resulting from the slat-cove 
separated flow are contributing factors. Certainly inade- 
quate turbulence modeling in the wake may be a contrib- 
uting factor as well, but there is no evidence at this time 
to point to its being a primary source of error. Determin- 
ing the actual causes will be difficult because there are 
almost no experimental flow field data over the slat it- 
self at this high Reynolds number. It will also be diffi- 
cult to draw conclusions about the ability of turbulence 
models to accurately predict processes such as slat wake 
/ boundary layer merging until we are confident that the 
slat flow field, and consequently its initial wake profile, 
are computed adequately to start with. 

6.2.2 Effect of Turbulence Model 

For the remainder of the computations in this paper 
(all at 19°), the ending transition locations from the ex- 
periment - as given in Table 1 - are used, with the ex- 
ception mentioned earlier that the slat upper surface 
transition is placed at x/c = -0.0847 to avoid laminar 
boundary layer separation. No boundary layer suction is 
employed. However, it is recognized that the flow over 
the slat is not being modeled correctly; therefore the re- 
sults are deficient in that regard. 

Also, all remaining computations are performed on 
grids with the tunnel walls modeled. The effect of the 
tunnel walls using EASM is shown in figures 16(a) - 
(d). The largest effect is on the profiles on the flap. 



q/a M 

(a) x/c = 0.1075 (main) 



q/ a oo 

(b) x/c = 0.45 (main) 


Fig. 15. Effect of slat transition location on 30P-30N velocity 
profiles, S-A model on free-air grid, alpha=8°. 


where including the walls causes the wakes to be shifted 
upward because the streamlines are constrained from 
turning as much. The wake deficits are also, for the most 
part, deeper when the walls are modeled. These wall ef- 
fects are consistent with those reported by Cao et al. 22 

Results using the S-A model and EASM are compared 
in figures 17(a) - (d) and in figures 18(a) - (d). The ve- 
locity profiles predicted by the two models are very sim- 
ilar, as are the turbulent shear stress profiles, which are 
rotated so as to be taken with respect to the directions 
parallel (u ) and normal ( v' ) to the body surface. How- 
ever, EASM tends to give slightly lower peak «V lev- 
els over the flap, in general. Also, at the first station on 
the flap, EASM predicts the peak w V to be deeper in 
the flap boundary layer than the S-A model. The reason 
for this shift is currently not known. In spite of this dif- 
ference, there is very little impact on the computed ve- 
locity profiles. 

Computations with the S-A model require roughly be- 
tween 3000 - 6000 multigrid cycles to converge, de- 
pending on the case and initial conditions, at a CPU cost 
of approximately 10 \is / grid point / multigrid cycle on 
a single-processor CRAY C90. EASM is more depen- 
dent on its initial conditions, and therefore can often re- 
quire more multigrid cycles to reach convergence, at a 
CPU cost of approximately 13 \is / grid point / multigrid 
cycle. 
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9/a- 

(a) x/c = 0.45 (main) 



(b) x/c = 0.89817 (flap) 



(c) x/c = 1.0321 (flap) 



9/a. 

(d) x/c = 1.1125 (flap) 




(b) x/c = 0.89817 (flap) 



(c) x/c = 1.0321 (flap) 



9/a. 

(d) x/c = 1.1125 (flap) 


Fig. 16. Effect of walls on 30P-30N velocity profiles, EASM, Fig. 17. Effect of turbulence model on 30P-30N velocity profiles, 

transition set, alpha=19°. transition set, walls modeled, alpha=19°. 
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uV 

(b)x/c = 0.89817 (flap) 



u'v' 

(c)x/c= 1.0321 (flap) 



uV 

(d)x/c= 1.1125 (flap) 


7. CONCLUSIONS 

The flows over two different multi-element airfoil 
configurations have been computed using linear eddy 
viscosity turbulence models and a nonlinear explicit al- 
gebraic stress model. A subset of recently-measured 
transition locations using hot film on a McDonnell Dou- 
glas configuration has been presented, and the effect of 
transition location on the computed solutions has been 
explored. Based on the computed results and compari- 
son with experimental data, the following conclusions 
are made: 

1. Specifying transition location correctly is crucial to 
the accurate computation of boundary layer velocity 
profiles. Poor boundary layer profile predictions lead to 
poor trailing-edge near-wake profile predictions and, 
therefore, poor predictions in the far wake region. Thus, 
inadequate wake predictions cannot be attributed solely 
to inadequate turbulence modeling, but must be attribut- 
ed in part to deficiencies in transition prediction on the 
generating element. Recently measured transition loca- 
tions are helpful in this regard. Using them for the main 
element improves the prediction of its boundary layer 
thickness, skin friction, and wake profile shape. 

2. The flow field over the slat is not well understood 
and is difficult to predict. Using the measured transition 
locations for the slat either results in slat separation or 
else results in only marginal improvement in slat wake 
profiles in comparison with experimental data. Possible 
reasons for this disagreement are: (i) the fact that CFD 
does not model the transition process itself, which may 
be more important on the slat than on the other elements, 
(ii) the hot film array used to measure the transition lo- 
cations may have some effect on the flow field and pos- 
sibly the transition locations themselves, and (iii) 
important slat physics (such as possible 3D effects, lam- 
inar separation bubbles, and shear-layer transition and/ 
or unsteadiness in the cove region) are not being mod- 
eled in the computations. Future CFD efforts in this area 
may be hampered by the fact that detailed experimental 
flow field measurements can be difficult to obtain in the 
slat region, particularly at high Reynolds numbers be- 
cause of the smaller scales. However, since the compu- 
tation of the slat flow field appears to represent a key 
roadblock to successful predictions of multi-element 
flows, some experimental efforts to document the flow 
physics in this region may be prudent. 

3. Although there are some variations, the nonlinear 
explicit algebraic stress turbulence model overall gives 
very similar results to linear eddy viscosity turbulence 
models. 
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Fig. 18. Effect of turbulence model on 30P-30N turbulent shear 
stress profiles, transition set, walls modeled, alpha=19°. 
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